Challenging the paradigm of singularity excision in gravitational collapse 
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A paradigm deeply rooted in modern numerical relativity calculations prescribes the removal of those regions 
of the computational domain where a physical singularity may develop. We here challenge this paradigm by 
performing three-dimensional simulations of the collapse of uniformly rotating stars to black holes without ex- 
cision. We show that this choice, combined with suitable gauge conditions and the use of minute numerical 
dissipation, improves dramatically the long-term stability of the evolutions. In turn, this allows for the calcula- 
tion of the waveforms well beyond what previously possible, providing information on the black-hole ringing 
and setting a new mark on the present knowledge of the gravitational-wave emission from the stellar collapse to 
a rotating black hole. 
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Numerical relativity simulations have recently recorded im- 
portant breakthroughs, which have allowed for long-term sta- 
ble and accurate calculations of curved and highly dynam- 
ical spacetimes, thus increasing their potential of providing 
physically relevant predictions for gravitational-wave astron- 
omy. Behind this rapid progress are various novel approaches, 
some of which involve the dropping of assumptions or tech- 
niques that were considered to be important or simply neces- 
sary. Firstly, the dismissal of the 3+1 ADM formulation of the 
field equations, which, after large-scale efforts 1 1], has shown 
not to have the stability properties necessary for long-term 
fully three-dimensional (3D) simulations. In lieu of the ADM 
equations, new formulations have been either reconsidered (as 
in the case of the conformal and traceless formulation of the 
ADM equations |2]) or investigated for the first time in 3D 
simulations (as in the case of the harmonic formulation of the 
Einstein equations 0). Both approaches have been shown to 
provide long-term stability on timescales sufficiently large to 
evolve accurately a large class of spacetimes, including black 
holes OllilS] and neutron stars 0,01 ■ Secondly, the abandon- 
ing of the use of numerical grids with uniform spacing, in lieu 
of which several codes now use mesh-refinement techniques 
(either fixed |8] or with adaptivity |4]). This conceptually 
simple but technologically challenging improvement allows 
to concentrate computational resources where the truncation 
error needs to be the smallest, while saving them where they 
may not be needed. In turn, non-uniform grids have allowed 
to place the outer boundaries of the computational domain at 
very large distances, thus reducing the influence of inaccu- 
rate outer-boundary conditions and making it possible to have 
the wave-zone within the domain and to extract there the pre- 
cious gravitational-wave signal |9]. Thirdly, for some years 
now, successful long-term 3D evolutions of black-hole vac- 
uum spacetimes have been possible only thanks to the use of 
excision techniques (see, e.g., [ 10] for a first example), that 
is by ignoring the spacetime regions inside black-hole hori- 
zons. These are causally disconnected from the outside and 
should have no effect on the rest of the evolution as long as 



a suitable treatment of the equations is made at the excision 
surface. Furthermore, recent simulations of the collapse of 
rotating neutron stars to Kerr black holes 0,(01 have shown 
the effectiveness of excision techniques also for spacetimes 
with matter, where they were applied separately to the field 
equations and to the hydrodynamical equations. In those sim- 
ulations, in fact, the use of excision has extended considerably 
the lifetime of the simulations, allowing for an accurate inves- 
tigation of the dynamics of the trapped surfaces formed dur- 
ing the collapse and for the extraction of the first gravitational 
waveforms from 3D collapse to rotating black holes 01 . 

Although the assumption that a region of spacetime that is 
causally disconnected can be ignored without this affecting 
the solution in the remaining portion of the spacetime is cer- 
tainly true for signals and perturbations travelling at physical 
speeds, numerical signals, such as gauge waves or constraint 
violations, may travel at velocities larger than that of light and 
thus leave the physically disconnected region. Indeed, this 
is what is commonly observed when excising a topologically 
spherical surface in Cartesian coordinates within a conformal 
and traceless formulation of the Einstein equations and with- 
out the massive use of numerical dissipation as in ref. 1 3]. 

In this paper we challenge the paradigm of singularity exci- 
sion and show that accurate numerical simulations can be car- 
ried out even in the absence of an excised region and indepen- 
dently of whether the spacetime is vacuum and the singularity 
modelled as a "puncture" 00]. In addition, we show that the 
absence of an excised region improves dramatically the long- 
term stability, allowing for the calculation of the gravitational 
waveforms well beyond what previously possible and past the 
black-hole quasi-normal-mode (QNM) ringing. 

The simplest and most impressive way of proving the im- 
provements resulting from not excising the region inside the 
apparent horizon (AH) is by reconsidering the same initial 
data of ref. |9], which lead to gravitational collapse to rotat- 
ing black holes. More specifically, we consider rotating rel- 
ativistic stars calculated as equilibrium solutions of the Ein- 
stein equations, with a polytropic equation of state p = Kp T , 
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with r = 2 and with the polytropic constant initially set to 
K lD = 100 1 23]. To illustrate the collapse for the range of 
neutron stars rotating either very slowly or near the mass- 
shedding limit, we consider two representative models, indi- 
cated as Dl and D4 in ref. lHUl . with different initial angu- 
lar momenta. The first one is a slowly rotating star of mass 
M = 1.67 Mq, circumferential equatorial radius R e = 11.43 
km and a = J/M 2 = 0.21, where J is the angular momen- 
tum; the second one has with M = 1.86 M©, R e = 14.25 km 
and is rotating close to the mass-shedding limit with a = 0.54. 

The numerical methods used for these evolutions are the 
same as discussed in ref. |9[. In particular, we use a con- 
formal and traceless formulation of the ADM equations 
which are solved with fourth-order finite-difference opera- 
tors, while we evolve the hydrodynamical equations with 
the Whisky code 1 1 311 - implementing high-resolution shock- 
capturing (HRSC) techniques with a variety of approximate 
Riemann solvers and reconstruction algorithms and an over- 
all second-order truncation error |11]. The numerical grid 
setup makes use of the same mesh refinement implemented 
in the numerical infrastructure described in ref. |&|]. Besides 
a fixed number of refinement levels which are present already 
on the initial slice, new refined levels are added at predefined 
positions during the evolution. More specifically, as the star 
collapses, the innermost (most refined) grid-box is progres- 
sively further refined with box-in-box grids, the final one hav- 
ing a resolution of Ax ~ 0.02 M ; the outermost grid, instead, 
which is not further refined, has a resolution of Ax ~ 1.5 M, 
sufficient to capture the details of the gravitational radiation. 
In this way, our outer boundaries are placed at ~ 165 M, so 
as to minimize the influence of our imperfect outer-boundary 
conditions on the very small gravitational-wave signal. We 
note that the "switching-on" of different levels of resolution 
does have a small but appreciable effect on its dynamics, es- 
sentially due to the necessary interpolation to the new refined 
grids, which also changes the growth-time of the most unsta- 
ble, exponentially growing mode. This effect is negligible in 
the evolution of the matter but can influence the waveforms. 
To provide a direct comparison with 0], we have here fol- 
lowed the same approach. 

Fig. ^ summarises the dynamics of the matter for model 
D4 by showing the time evolution of the maximum of the 
rest-mass density /5 max (continuous lines) and of the total rest- 
mass (dashed lines) when normalized to their initial val- 
ues (a similar behaviour is observed also for the slowly ro- 
tating model Dl). The evolution obtained without excision is 
to be contrasted with the one in which the excision is made 
(the used excision technique is described in detail in ref. Q) 
and which is indicated with thick lines. In that case, the ex- 
cision was started soon after an AH fl5ll was found (this is 
marked with circles on the two curves) and the correspond- 
ing evolution terminated at t > 77 M, when the code crashed 
with large violations in the Hamiltonian constraint. Note that 
while the maximum values for the rest-mass density attained 
during the evolution are comparable in the two cases, the evo- 
lution without excision can be carried out to much later times 
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FIG. 1: Evolution of the normalized maximum of the rest-mass den- 
sity Pmax (continuous line) and of the total rest-mass M, (dashed 
line) for model D4. Circles indicate the times at which the AH is 
first found and thick lines the results obtained with excision. 

(i.e. t > 300 M ) without appreciable loss of accuracy or sign 
of instability. Indeed, as the matter collapses and concentrates 
over a very few gridpoints (and ultimately on only one), the 
high accuracy of the HRSC methods is able to conserve the 
rest mass to very high precision with a loss of less than 0.2% 
up to when the rest-mass density distribution is diffused as a 
result of the poorly resolved gradients. 

While very little extra is needed for the evolution of the 
hydrodynamical quantities in the absence of an excision algo- 
rithm, a stable evolution of the Einstein equations requires at 
least two important ingredients. The first one is represented 
by gauge conditions for the lapse function a with suitable 
singularity-avoidance properties. Our experience has shown 
that hyperbolic X-driver slicing conditions of the form dtot = 
—/(a) a 2 (K — Kq) , with / > and K being the trace 
of the extrinsic curvature at t = 0, are essential to "freeze" 
the evolution in those regions of the computational domain in- 
side the AH, where the metric functions experience the growth 
of very large gradients. In practice, we confirm that using 
the generalized "1+log" slicing condition 11611 as obtained by 
setting / = 2/ a provides the desired singularity avoidance 
and is computationally efficient. With a good choice for the 
slicing condition, the results do not depend sensitively on the 
gauge condition for the shift. We have found that the use of 
the "Gamma-driver" shift conditions discussed in ref. II ill is 
sufficient to compensate for the "slice-stretching" induced by 
the singularity-avoiding slicing. 

The second important ingredient is the introduction of an 
artificial dissipation of the Rreiss-Oliger type 1 17] on the 
right-hand-sides of the evolution equations for the spacetime 
variables and the gauge quantities (No dissipation is intro- 
duced for the hydrodynamical variables.). This is needed 
mostly because all the field variables develop very steep gra- 
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FIG. 2: Evolution of the L2 norm of the Hamiltonian-constraint vio- 
lation for model Dl. Different curves refer to the violation computed 
only inside the AH (short-dashed line), only outside it (continuous 
line) or when the excision is made (long-dashed line in the inset). 
The circles show the rescaled violation for a higher resolution. 

dients in the region inside the AH. Under these conditions, 
small high-frequency oscillations (either produced by finite- 
differencing errors or by small reflections across the refine- 
ment boundaries) can easily be amplified, leave the region 
inside the AH and rapidly destroy the solution. In practice, 
for any time-evolved quantity u, the right-hand-side of the 
corresponding evolution equation is modified with the intro- 
duction of a term of the type £ diss (u) = —eAx^d^.u, where 
e is the dissipation coefficient, which is allowed to vary in 
space. We have experimented with configurations in which 
the coefficient was either constant over the whole domain or 
larger for the gridpoints inside the AH. No significant differ- 
ence has been noticed between the two cases. Much more 
sensitive is instead the choice of the value of e. In the sim- 
ulations reported here we have set e = 0.01 for model Dl 
and £ = 0.0075 for D4, respectively. However, smaller val- 
ues (e.g. e — 0.005) are not sufficient to yield the long-term 
stability discussed here, while larger values (e.g. e = 0.05 
and e = 0.01, respectively) alter significantly the waveforms, 
which would not match the ones obtained without dissipation 
up until the the latter can be computed. 

In the region within the AH, as a consequence of large gra- 
dients in the spacetime variables, which cannot be resolved 
sufficiently despite the use of several mesh-refinement levels, 
the solution of the Einstein equations becomes increasingly 
inaccurate as the collapse proceeds. Fig.|2]offers a measure of 
this loss of accuracy by showing the time evolution of the L2 
norm of the Hamiltonian-constraint violation for model Dl. 
To distinguish the different amplitudes of the errors, the vio- 
lation has been computed separately for the domain inside the 
AH and for the rest of the grid. 

Several aspects of Fig. [2] are worth noting. First, the two 



errors are considerably different in amplitude, the one inside 
the AH being much larger. Second, the sudden decrease in 
the long-dashed and continuous lines at the time of AH for- 
mation at t ~ 60 M is due to the fact that some gridpoints 
(those with the largest violations) are no longer included in 
the calculation of the L2 norm. Third, the violation outside 
the AH grows only at a slow pace over the remaining evolu- 
tion. Finally, shown with circles is the violation for a high- 
resolution simulation. The values are rescaled to second or- 
der before AH formation (when the dominant truncation er- 
ror comes from hydrodynamics and the HRSC methods used 
are only second-order) and to third order afterwards (when 
the smooth accretion flow of the atmosphere boosts the accu- 
racy of the hydrodynamics to third-order). No convergence 
of the violation is seen inside the AH. Overall, Fig. |2] shows 
that the solution of the Einstein equations inside the AH is by 
and large incorrect, but these errors remain confined within 
the AH and do not contaminate the overall accuracy of the 
simulation. Shown for comparison in the inset is the violation 
computed when the singularity is excised, showing an expo- 
nential growth soon after the AH is found. 

Note that the ability to perform these long-term evolutions 
cannot be related in any manner to the use of a "puncture" pre- 
scription for the physical singularity, as recently done in simu- 
lations of binary black holes |0,|5|. However, it is possible that 
the stability provided here by the singularity-avoiding gauge 
and made numerically more robust by the use of dissipation is 
closely related to the one discussed for punctures 1 18]. 

Besides a long-term stability and the possibility of fol- 
lowing the collapse well beyond what was possible with 
the use of excision techniques, the most dramatic advan- 
tage produced by the approach suggested here can be ap- 
preciated in the study of the gravitational radiation produced 
during the collapse. As in ref. we have extracted the 
gravitational-wave information through an approach in which 
the spacetime is matched with the non-spherical perturbations 
of a Schwarzschild black hole described in terms of gauge- 
invariant odd Q)°* and even-parity metric perturbations, 
where £, m are the indices of the angular decomposition. In 
Fig. 13 we report the lowest-order multipoles Q^ = -V^o, 
where A = y/2(l + 2)\/(l-2)\. The left panel, in partic- 
ular, shows the signal detected by an observer at a distance 
of 42 M for the slowly rotating model Dl; the right panel, 
instead, shows the same multipole extracted at 37 M for the 
rapidly rotating model D4. In both cases the thick lines re- 
fer to the evolutions carried out with the excision technique 
and terminate when the code crashed (cf. left panel of Fig. 1 
in ref. [9]). This comparison shows that it is now possible 
to detect the complete gravitational-wave signal produced by 
the collapse of a relativistic star to a rotating black hole well 
beyond what was previously possible in either 2D [ 19] or 3D 
simulations |9]. The estimated error on the phase is < 1% and 
< 5% on the L2 norm of the amplitude. 

Fig 13 also highlights that the signal can be rather different 
both in amplitude and form in the two cases with the excep- 
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FIG. 3: Lowest gauge-invariant multipole for a slowly rotating star (left panel) and a rapidly rotating one (right panel). Thick lines refer to 
evolutions carried out with excision and terminate at a code crash, while the circles indicate the time of the AH formation time. 



tion of the final parts, when the signal is dominated by the 
black-hole QNM ringing. This happens between ~ 80 M and 
~ 120 M, where the signal we extract matches very well with 
a perturbative one computed using the frequencies given in 
ref. HJ for the angular momenta and masses of our models. 
The very good agreement with the perturbative results is an 
important confirmation of the accuracy of our results, which 
are intrinsically plagued by the extreme weakness of the emit- 
ted gravitational radiation. Furthermore, the richness of de- 
tails in the two waveforms opens the prospects that a careful 
characterization of the waveforms will provide important in- 
formation on the properties of the black hole as well as on 
those of the collapsing matter |21]. 

A straightforward analysis of the now-complete 
gravitational-wave signal computed in these simulations 
allows to improve the estimates provided in ref. |9] for the de- 
tectability of the gravitational collapse of a uniformly rotating 
poly tropic star at a distance of 10 kpc. More specifically, the 
energy efficiency in the emission of gravitational radiation 
is E D1 /M = 3.3 x 10~ 7 and E V JM = 3.7 x 10~ 6 re- 
spectively, with an overall accuracy of ~ 10%. The resulting 
signal-to-noise ratios are then: (S/N)^^ ~ 0.27 - 2.1, 

(S/ivCT c 1.2 - 11, and (S/N)lT_ Di c 3.3 - 28for 
detectors such as Virgo/LIGO, advanced LIGO or Dual |22|. 

In conclusion, we have presented the first 3D calculations of 
the gravitational collapse of uniformly rotating stars to black 
holes without excision. This choice, together with suitable 
gauge conditions and the use of minute numerical dissipa- 
tion, improves dramatically the long-term stability of the evo- 
lutions, providing the most accurate waveforms of this pro- 
cess to date. While our approach represents a challenge to the 
paradigm of singularity excision, it does not necessarily imply 
that all excision techniques should be expected to lead to in- 



stabilities. Rather, it highlights that, for a conformal traceless 
formulation of the Einstein equations and in highly dynami- 
cal spacetimes, the excision of a spherical surface in Cartesian 
coordinates may be more of a problem than a solution. 
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